### to replicate Table 2 of Cao and Ward: Transnational Climate Governance Networks and Domestic Regulatory Action; 
### for a special issue in International Interactions;
### Thursday, 2.June.2016
### XC at Penn State: email xuc11@psu.edu for questions. 

dd <- c("C:/XunCao/Research/TCG/final/data code/table 3")  ### you need to set up your own working directory though. 
setwd(dd)

set.seed(300)

options(scipen=200)
library(xtable)

### 1: run the model using gbme.r from Peter Hoff: 
Yd<-dget("dat.tcg3")[[1]] ## this is the matrix for the dependent variable which counts the number of shared TCGs between countries;
Y<-Yd
n<-dim(Yd)[1]             ## you need to specify n to run the gbme.r
Xd<-dget("dat.tcg3")[[2]] ## this is N by N by 16 an array for dyadic covariates; 


source("gbme.r")          ## source the code for gbme.r;
gbme(Y=Y,Xd=Xd,fam="poisson",k=3,odens=100, directed=F, NS=200000, oround=6)       # run; this will generate two output files ==> this is going to take a while ... maybe more than one hour;
                                                                                   # you need OUT for the posterior distributions; 


    
      
### 2: summarize the results: after running the MCMCs;
out<-read.table("OUT",header=T)             #read in output
burn<-500 # the length of burn in: 500 maybe is a little too much ... but we have 2000 rows in OUT anyways; 

tit<-c("trade proxy", "distance", 
        "domestic legislation", 
"fossil",
         "green NGOs",
#         trad[country.set,country.set],
         "trade_share", 
         "real gdp pc", 
#         "polity", 
"civil liberty",
"veto player",
         "cci", 
#         oilgascoal, 
         "CO2 pc", 
         "EU", 
          "bank deposits per gdp", 
         "federation", "population", "iso14000 per gdp" #, 
         #"common language", "polity similarity"
         )
       

# we need some summary statistics of the posteriors:
Sum.post<-NULL
for (j in c(20, 4:19, 21:dim(out)[2], 3)){
                               sum.post<-c(quantile(out[,j][-c(1:burn)], c(.025,.05, .95, .975))[1], quantile(out[,j][-c(1:burn)], c(.025,.05, .95, .975))[2],
                               mean(out[,j][-c(1:burn)]), 
                               quantile(out[,j][-c(1:burn)], c(.025,.05, .95, .975))[3],quantile(out[,j][-c(1:burn)], c(.025,.05, .95, .975))[4], 
                               sd(out[,j][-c(1:burn)]))
                               Sum.post<-rbind(Sum.post, sum.post)
                               }

Sum.post<-round(Sum.post, digits=4)
colnames(Sum.post)<-c("2.5%", "5.0%", "Mean", "95.0%", "97.5%", "St.Dev")
rownames(Sum.post)<-c("constant", tit, colnames(out)[c(21:dim(out)[2], 3)])
Sum.post
xtable(Sum.post, digits=c(rep(4,7)))
